function dM=Blocheqs(t,M,B,theta,omega)
val=parameter();
Rp=val.Rp;
T1=val.T1;
T2=val.T2;
s=val.s;
%Me0=norm(s)*Rp/(Rp+1/T1);
 
%   % square   
  
  dM=gamma_e*cross(B,M)+...
     [-(cos(theta)^2*M(1)-cos(theta)*sin(theta)*M(3))/T2;...
     -M(2)/T2;...
     -(sin(theta)^2*M(3)-cos(theta)*sin(theta)*M(1))/T2]+...
     -M/T1+Rp*square(omega*t)*(s-M); 
end